Comprehensive study of the critical behavior in the diluted antiferromagnet in a field 
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We study the critical behavior of the Diluted Antiferromagnet in a Field with the Fethered Monte 
Carlo formalism. We compute the critical exponents (including the elusive hyperscaling violations 
exponent 6). Our results provide a comprehensive description of the phase transition and clarify the 
inconsistencies between previous experimental and theoretical work. To do so, our method addresses 
the usual problems of numerical work (large tunneling barriers and self-averaging violations). 
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Understanding collective behavior in the presence of 
quenched disorder has long been one of the most chal- 
lenging and interesting problems in statistical mechanics. 
One of its simplest representatives is the random field 
Ising model (RFIM), which has been extensively stud- 
ied both theoretically and experimentally^ The RFIM 
is physically realized by a diluted antiferromagnet in an 
applied magnetic field (DAFF). 

It is known that the D = 3 DAFF/RFIM undergoes 
a phase transition, but the details remain controversial, 
with severe inconsistencies between analytical, experi- 
mental and numerical work. A scaling theory is gen- 
erally accepted, where the dimension D of the system is 
replaced by D — 9 in the hyperscaling relation. This third 
independent critical exponent, believed to be ~ 1.5, is 
inaccessible both to a direct experimental measurement 
and to traditional Monte Carlo methods. 

The values of the remaining critical exponents, seem- 
ingly more straightforward, are also controversial. On the 
experimental front, different ansatze for the scattering 
line shape yield mutually incompatible estimates of the 
thermal critical exponent, namely v = 0.87(7) (Ref. 0), 
or v = 1.20(5) (Ref. i). ' Furthermore, the experimen- 
tal estimate of the anomalous dimension, rj = 0.16(6) 
(Ref. d), violates hyperscaling bounds, if one is to be- 
lieve the experimental claims of a diverging specific heat 
(a>0)4 

On the other hand, the numerical determination of 
v has steadily shifted, the most precise estimate being 
1.37(9) (Ref. [g), inconsistent with the experimental val- 
ues and barely compatible with a ~ 0. The value of a 
itself is very hard to measure in a numerical simulation. 6 

More fundamentally, the smallness of the magnetic ex- 
ponent P, combined with the numerical observation of 
metastability^ has led some authors to suggest that the 
transition in the DAFF may be of first order. 

Ultimately, the physical reasons for this confusion be- 
tray the fact that the traditional tools of statistical me- 
chanics are ill-suited to systems with rugged free-energy 
landscapes. Both experimentally and numerically, the 
system gets trapped in local minima, with escape times 
that grow as log r ~ £ e (£ is the correlation length). This 
not only makes it exceedingly hard to thermalize the sys- 
tem, but also generates a rare-events statistics, causing 



self-averaging violations^ 

In this letter we study the DAFF with the Tethered 
Monte Carlo (TMC) formalism,- Our approach restores 
self-averaging and is able to negotiate the free-energy bar- 
riers of the DAFF to equilibrate large systems safely. It 
also provides direct access to the key parameter 9. We 
thus obtain a comprehensive picture of the phase transi- 
tion, consistent both with analytical results for the RFIM 
and with experiments on the DAFF, and shed light on 
the reasons behind the previous discrepancies. 

In the following we provide a brief outline of the teth- 
ered formalism applied to the DAFF (see Refs. and 
[Tol for details). We note, however, that we give most of 
our physical results translated into the familiar canon- 
ical language. In a tethered computation, we run sim- 
ulations where one (or more) order parameters of the 
system are (almost) constrained. In this way, we elimi- 
nate the need for exponentially slow tunneling caused by 
the free-energy barriers associated to these parameters. 
From these tethered simulations the Helmholtz effective 
potential is accurately reconstructed with a fluctuation- 
dissipation formalism. 

We consider a system with N — L D spins, s x = ±1, 
on the nodes of a cubic lattice with periodic boundary 
conditions and interacting through the Hamiltonian 



H = 



-hM—h s M B = U-hM-h s M s . (1) 



Here h and h s are the applied fields, coupled to the mag- 
netization and staggered magnetization, 

M = Nm = J2^ s ^ M.= ^ e^e 1 * (2) 



We are ultimately interested in h s = 0, but we will find 
this parameter useful. The quenched occupation vari- 
ables e x are 1 with probability p — 0.7 and zero other- 
wise (this value is chosen to be far both from the perco- 
lation threshold and from the pure system). For D = 3, 
the system undergoes a paramagnetic-antiferromagnetic 
phase transition, where m s is the order parameter. 

Let us consider a single sample of the system (i.e., 
a fixed {e^})- In our tethered computation, we define 
smooth magnetizations rh and rh s by coupling m and m s 



2 



to Gaussian baths and work in a statistical ensemble for 
fixed (m, m 8 ) with weight^ 

u>(m, m s ; {s x }) oc e' 0u j(m, m)~f(m a , m a ), (3) 

where j(x,x) = e N ^ x ~ x \x — x)( N ~ 2 ^ 2 0{x — x), and 
0(x — x) is the step function. The smoothing procedure 
shifts the mean value of the parameters, sox~a; + l/2. 
This ensemble is related to the canonical one through 
a Legendre transformation. For instance, the partition 
function of the system is 

Z = J dmdm s w(m> s ; {s x }) e l3N< - hAl+h ^ 

M (4) 

-/**.•-""•'"•—" 

where i?Ar(m,7h s ) is the Helmholtz effective potential. 

We can reconstruct 17jv from computations at fixed 
(rh,rh s ) via the so-called tethered field (b,b s ) 

s=1 _l/2-l/iV k = 1 _lj2-W 
m — m m s — m s 

In particular, the gradient VJ?/v is 

(df2 N /dm, dQ N /dm s ) = ((6)m,™ s , (b s )rh,rh a ) ■ (6) 

The notation (• • • )m,m a denotes tethered expectation val- 
ues, computed with weight ([3]). 

A TMC computation consists in a set of independent 
Monte Carlo simulations at fixed (m, m s ) that are then 
combined to reconstruct Qn- Note that the effective po- 
tential (as a function of the magnetizations) has all the 
information about the system in the tethered ensemble, 
just as the free energy (as a function of the applied fields) 
has all the information in the canonical ensemble. 

The canonical averages at fixed (h, h s ) can be recovered 
with Eq. (UJ. Note that, according to ©, this integral is 
dominated by saddle points (m, m s ) such that 

(b)rh,rh s = 0h, (b a }m,rh, = PK. (7) 

We can determine the relative weights of different sad- 
dle points by line-integrating the tethered field along any 
connecting path. We are interested in the case h s = 0. 

So far we have summarized the application of TMC for 
a single sample. Since it consists of simulations at fixed 
(m,m s ), it eliminates the need to tunnel between coex- 
isting phases and, hence, equilibrates the system much 
faster than a canonical simulation. However, we still face 
the serious problem of self-averaging violations. In prin- 
ciple, the definition of quenched disorder implies recon- 
structing the free energy with Q before computing the 
disorder average. In this work, however, we sample aver- 
age the Helmholtz potential rather than the free energy 
(a similar approach was taken in Ref. fill ). 

In order to motivate this approach, let us consider Fig- 
ure Q] — Top. We compare the tethered average (6 s )m,m 3 
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FIG. 1. (color online) Top: Tethered field (S S ) A ,™ 3 , Eq. ©, 
at T — 1.6 and m = 0.11, for two individual samples of an 
L = 24 system (□ and ■) and for the sample average (•) 
as a function of m s . The field is self-averaging in the region 
outside the two external zeros. The errors cannot be seen 
at this scale. Bottom: Effective potential J7jv(m = 0.11, m s ) 
obtained by integrating the averaged tethered field of the top 
panel. The two antiferromagnetic minima are separated by 
a very large barrier (the escape time is r ~ exp[NAO]), and 
there is no paramagnetic minimum. 

for two individual samples with the disorder average over 
1000 samples. The zeros of this latter curve separate an 
internal gap with chaotic fluctuations, where the field 
vanishes in the thermodynamical limit, from an external 
region where the field is actually self-averaging. 

We exploit the situation by considering a small, but 
finite, value of h s . The saddle point defined by this field 
will be in the self-averaging region. We can therefore 
solve the saddle-point equations on average, rather 
than sample by sample. Only afterwards do we make 
h s — > in the solution (this is analogous to the math- 
ematical definition of spontaneous symmetry breaking). 
The limit h s — + is essentially equivalent to considering 
a 'smeared' saddle point and averaging over all rh s 

70) fn = J dm s M AA e-*( A '^- fi °] . (8) 

S7q is a normalization constant. Since we work at fixed 
m, On is just the one-dimensional integral of (bs)^ A . 

The other saddle-point equation, (b) rh — /3h, defines a 
one-to-one relation rh(h) so that (O)^^ and the canoni- 
cal (0)(h) both tend to the same thermodynamical limit 
(ensemble equivalence). Furthermore, for finite lattices 
{0) 1% is better behaved statistically and arguably more 
faithful to the physics of an experimental sample. There- 
fore, we shall identify (0)(h) = and use the more 
familiar canonical notation. See Refs. HI and[I3 for a more 
detailed study of this ensemble equivalence. 

We have used the above outlined procedure to ther- 
malize the DAFF for temperatures down to T = 1.6 and 
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FIG. 2. (color online) Top: Correlation length £/L as a func- 
tion of the applied magnetic field h for T = 1.6. The curves 
intersect, marking a second-order phase transition. Bottom: 
Scaling plot of £ as a function of T for h = —2.13, showing 
large corrections to leading scaling (we use v — 1.05). 

sizes up to L = 32 (1000 samples for L = 8, 12, 16, 24 and 
700 samples for L = 32). For each size we simulate a grid 
of ps 150 points in the (m,m s ) plane (5 values of m, and 
fa 30 values of m s on each). We also use temperature 
parallel tempering. This is only necessary to thermalize 
L > 24, but it is convenient for smaller lattices because 
we are also interested in the T dependence. Thermaliza- 
tion is ensured using the methods described in Ref. Il2l . 
We provide more technical details in Ref. O 

The first interesting physical result is the effective po- 
tential itself. Some authors have found metastable be- 
havior in the DAFF, interpreted as a sign of a first-order 
transition.- This should manifest as the coexistence of 
antiferromagnetic and paramagnetic minima in fl. How- 
ever, see Figure Q] — Bottom, our results exhibit only two 
antiferromagnetic minima, separated by a very large free- 
energy barrier. In a canonical simulation, the system 
tunnels back and forth between the two, with an escape 
time r ~ exp[AAi?]. This explains the metastable be- 
havior observed in previous work (and the difficulty to 
thermalize large samples with canonical methods) , but is 
inconsistent with a first-order scenario. 

Of course, we could be looking at a value of rh (equiv- 
alently, of h) far from the critical point. In order to 
find the phase transition, we compute the usual second- 
moment cor relation lengt h We use the propagator 
Fh(k) = N(<fi(k)<f>(—k))(h), where <f> is the staggered 
Fourier transform of the spin field. 

We have plotted £(h)/L at T = 1.6 as a function of 
the applied field h in Figure [5] — Top. The curves for dif- 
ferent L show very clear intersections, marking the onset 
of a second-order phase transition. In order to estimate 
the critical exponents, we apply the quotients method^ 
We consider the ratios of physical observables for system 



L h*(L) Pjvh ot/vh vt 

8 -2.178(4) 0.0125(7) 0.887(5) 0.0765(25) 1.07(9) 

12 -2.140(5) 0.0104(5) 0.790(9) 0.0781(27) 1.01(4) 

16 -2.123(3) 0.0119(4) 0.742(7) 0.224(4) 1.10(15) 



TABLE I. Computation of the critical exponents using the 
quotients method. We extract our estimates from ratios of 
physical observables for sizes (L,2L), computed at the inter- 
section point of £/L. The first four columns give results for 
fixed T — 1.6 and the last one at fixed h — —2.13. 



L 


AF/N 


Fit range 


e 


X 2 /d-o.f. 


8 


0.03382(29) 


L > 8 


1.448(9) 


5.56/3 


12 


0.01756(15) 


L > 12 


1.469(13) 


0.44/2 


16 


0.01138(9) 


L > 16 


1.461(20) 


0.16/1 


24 


0.00608(5) 








32 


0.00392(5) 









TABLE II. Computation of the hyperscaling violations expo- 
nent 6 from the free-energy barriers AF. We report fits to 
AF — AL , for different ranges, giving the x 2 and the de- 
grees of freedom of each fit. Our preferred final estimate is 
9 = 1.469(20), taking the central value of the fit for L > 12 
and the more conservative error of the fit for L > 16. 



sizes (L,2L) } computed at the intersection point h*(L) 
of their respective £(h)/L. We have applied this method 
to d h i ~ L x +V^ and (mf)(/i) - i 2 * 9 /^" 3 in Table H 
Note that our estimate for (3 is very low, in accordance 
with previous numerical and experimental work. 

We can also estimate v from the temperature depen- 
dence of £ at fixed h, obtaining a second estimate ut 
(Table |T|. Both determinations of v should coincide, but 
we obtain Vh ~ 0.75 and vt ~ 1.05. We can see in Fig- 
ure [2] — Bottom that this discrepancy is due to strong 
scaling corrections. If one attempts a collapse of the 
curves, focusing on different ranges for £/L, the corre- 
sponding values of v vary from v ps 0.75 to v > 2, which 
explains the wide range of variation in previous numer- 
ical estimates of v. By safely locating the critical point 
and using the quotients method, we have minimized the 
scaling corrections, but not eliminated them completely. 

We need an additional critical exponent in order fully 
to characterize the critical behavior of the DAFF. This is 
the hyperscaling violations exponent 0, which can be re- 
lated to the free-energy barrier between the ordered and 
the disordered phase: AF oc L e } 4 - The computation of 
these barriers is very difficult with traditional methods, 
but straightforward with TMC. Indeed, we can identify 
AF with the A/2jv between the two saddle points (disor- 
dered and antiferromagnetic) defined by the critical h c . 

We can compute this barrier simply by evaluating the 

line integral of «S) A)As - /3h c , (M A A J along a path 
joining the two saddle points. We know that one of them 
will lie on the line m s = 0.5 (m s ~ 0). Therefore, we 
first integrate from the antiferromagnetic saddle point to 
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ih a = 0.5 at fixed m. We then integrate at fixed rh s = 
0.5 until we reach the disordered saddle point. We give 
the resulting values of Af2 N = AF/N in Table HH Our 
final estimate is 8 — i. 469 (20), incompatible with the 
8 = D — 1 of a first-order phase transition. 

Notice that the hyperscaling relation 2 — a = v(D — 8), 
coupled with our values for v and 8, predicts not only 
a divergence of the specific heat, as observed in experi- 
ments, but also a positive a . We could test this result di- 
rectly by computing C = (m) . Unfortunately, the quo- 
tients method is ill-suited to this quantity, whose scaling 
is more aptly described as C ~ A + BL a l v J£ Therefore, 
one needs extremely large values of L to reach the asymp- 
totic regime C ~ L a ' v . The behavior of the quotients in 
Table H] is consistent with this expectation. 

It has been proposed that 8 is not independent, but 
given by 8 = D/2 - /3/vM Combining Tables U and HI 
we see that our numerical results are indeed compatible 
with this two-exponent scenario. 

We can use our results to comment on the experi- 
mental situation. In an experimental study, the critical 
exponents are computed from fits to the scattering line 
shape S(k) — Sd(k) + S c (k), where the two terms distin- 
guish connected and disconnected contributions. In the 
two-exponent scenario, strongly supported by our data, 
the most singular term in Sd is the square of S c - This 
Ansatz was applied in Ref. |2j, yielding v — 0.87(7) and 
r\ = 0.16(6). Since rj = 8 — 1 + 2/3 /V, however, this last 
value violates hyperscaling bounds and is also incompat- 
ible with our results. Perhaps taking Sd = (S c ) 2 for the 
whole function, not just its singularity, is an excessive 
simplification. Clearly a better theoretical determina- 
tion of S(k) is needed. Our methods are well suited to a 



direct numerical approach to this question. 

We have used the tethered formalism to obtain a com- 
prehensive picture of the critical behavior of the DAFF, 
resolving the inconsistencies in previous work. This 
method restores self-averaging to the problem and is 
capable of handling rugged free-energy landscapes to 
equilibrate much larger systems than canonical paral- 
lel tempering. Our simulations show clear signs of a 
second-order phase transition and are consistent both 
with experiments on the DAFF and with analytical re- 
sults for the RFIM. The critical exponents 8 and f3/v 
(equivalently, r] and fj) are computed with a high pre- 
cision, although our simulations were not optimized for 
the computation of v (equivalently, of a). We obtain 
v = 0.90(15), consistent with a positive a. 

The tethered approach demonstrated in this paper has 
a very broad scope and we believe it can be fruitfully ap- 
plied to many systems featuring large free-energy barri- 
ers. Indeed, it has already been successfully implemented 
for hard-spheres crystallization^ Other promising av- 
enues are the study of Goldstone bosons and the equation 
of state for the D = 3 spin glass/-® or equilibrium and 
aging relaxation in a metastable phase (e.g., to prevent 
crystallization of supercooled liquids, see Ref. Ii~9h . 
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